依然滿滿的內容 先謝謝閱讀了XD
Day 22 和 Day 23，我們完成了 Crossing 這條線。
現在只要給一份 layout，evaluator 就能告訴我們：
每條 edge 的 crossing count
K
n_K
Phi
C
但 evaluator 只會打分數，不會產生答案。
真正開始跑 SA 之前，我們仍然需要替每個 node 準備第一組座標：
node 0 -> (x0, y0)
node 1 -> (x1, y1)
node 2 -> (x2, y2)
...
最簡單的方法，是在畫布上隨機撒點。
這個方法便宜、多樣，而且只要處理重複座標，就能很快產生一批 seeds。但它完全不看 graph topology：明明有 edge 相連的兩個 nodes，可能被放在畫布兩端；完全不相干的一群 nodes，反而擠在同一個角落。
graph topology
告訴我們誰和誰相連
random layout
完全沒有使用這份資訊
於是 SA 一開始不只要降低 crossings，還得先花很多步把 graph 的基本結構整理出來。
今天要介紹 Layout 的第一個版本：force-directed layout。
Force-directed layout 的直覺很像一套假的物理系統：
所有 nodes 彼此排斥
-> 避免全部擠成一團
有 edge 相連的 nodes 彼此吸引
-> 避免相連 nodes 離得太遠
每一輪計算合力並移動
-> 逐漸形成比較均勻的 graph shape
假設有一條 path：
0 -- 1 -- 2 -- 3 -- 4
隨機位置可能長成：
3 ●
0 ● ● 1
4 ●
● 2
Force-directed iterations 會嘗試讓相連的 nodes 靠近，同時讓所有 nodes 保持一些距離：
0 ●-----● 1-----● 2-----● 3-----● 4
這只是一個直覺圖。Force-directed layout 不是在模擬真實材料，也不知道什麼是 LCN；它只是用一組人為設計的 forces，把 topology 轉成座標。
Force-directed 是一整個演算法家族，不是唯一一條公式。
專案裡同時存在：
FR-style force model
公式較直接
適合當教學與 CPU fallback
OGDF-aligned FMMM force model
attraction 使用另一套 log2(d/L) 公式
還包含 multilevel、coarsening 與更完整的 iteration policy
Day 24 先使用 Fruchterman–Reingold，也就是 FR baseline。因為它最容易看出 repulsion、attraction 與 temperature 各自做什麼。
Day 25 再進入實際 FMMM GPU backend 的 tiling、shared memory、kernel fusion 與大圖路徑。
先把 FR 使用的理想距離記成 k：
area = width * height
k = k_scale * sqrt(area / N)
N 是 node 數量。畫布越大，理想距離越大；nodes 越多，每個 node 能分到的空間越小。
k_scale 則讓我們調整 layout 想要多鬆或多緊：
k_scale larger
-> nodes 傾向分得更開
k_scale smaller
-> graph 傾向更緊密
對 node i 和 node j：
dx = x[i] - x[j]
dy = y[i] - y[j]
distance = sqrt(dx * dx + dy * dy)
排斥力大小使用：
repulsive_magnitude = k * k / distance
方向則從 j 指向 i：
unit_x = dx / distance
unit_y = dy / distance
force_x = unit_x * repulsive_magnitude
force_y = unit_y * repulsive_magnitude
把大小與方向合併，也可以寫成：
force_x = dx * k * k / distance_squared
force_y = dy * k * k / distance_squared
距離越近，排斥越強；距離越遠，影響越小。
對一個 node i，需要累加其他所有 nodes 的貢獻：
for i in range(num_nodes):
force_x[i] = 0.0
force_y[i] = 0.0
for j in range(num_nodes):
if i == j:
continue
dx = pos_x[i] - pos_x[j]
dy = pos_y[i] - pos_y[j]
dist_sq = dx * dx + dy * dy
if dist_sq > epsilon:
force_x[i] += dx * k * k / dist_sq
force_y[i] += dy * k * k / dist_sq
這是一個標準的 all-pairs 計算：
node 0 -> read positions of 0 ... N-1
node 1 -> read positions of 0 ... N-1
node 2 -> read positions of 0 ... N-1
...
工作量是：
N * N
這一段通常會成為 force iteration 的主要成本。
如果兩個 nodes 剛好在相同位置：
dx = 0
dy = 0
distance = 0
直接帶進公式會除以零。因此實作通常會設定 epsilon：
if dist_sq <= epsilon:
continue
但這個做法有一個反直覺的副作用：兩個完全重疊的 nodes 得到的方向也是 (0, 0)，所以 repulsion 反而無法把它們分開。
更完整的實作可以對 coincidence 加入可重現的小擾動，或在後面的 quantization / repair 階段解決碰撞。不能只說「加 epsilon 就完成了」，因為數值安全與幾何合法性是兩個不同問題。
Repulsion 對所有 node pairs 生效；attraction 只沿著 edges 計算。
對 edge (u, v)：
dx = x[v] - x[u]
dy = y[v] - y[u]
distance = sqrt(dx * dx + dy * dy)
FR attraction 的大小是：
attractive_magnitude = distance * distance / k
距離越遠，拉回來的力量越強。
方向從 u 指向 v，所以同一條 edge 會對兩端產生相反的 force：
force on u += attraction toward v
force on v += attraction toward u
Naive CPU 寫法可以是：
for u, v in edges:
dx = pos_x[v] - pos_x[u]
dy = pos_y[v] - pos_y[u]
distance = sqrt(dx * dx + dy * dy)
if distance > epsilon:
magnitude = distance * distance / k
fx = dx / distance * magnitude
fy = dy / distance * magnitude
force_x[u] += fx
force_y[u] += fy
force_x[v] -= fx
force_y[v] -= fy
每條 edge 只處理一次，因此 attraction 的工作量是：
E
這裡已經出現兩種完全不同的資料形狀：
repulsion
dense all-pairs
cost is proportional to N * N
attraction
sparse edge list
cost is proportional to E
雖然它們最後都寫進 force_x[node]、force_y[node]，適合的 GPU mapping 並不相同。
把 repulsion 和 attraction 相加後，我們會得到每個 node 的 displacement direction：
disp[i] = repulsion[i] + attraction[i]
如果直接使用完整 force magnitude 更新位置，某些 nodes 可能一步衝過平衡點，下一輪又被拉回來：
left ---- equilibrium ---- right
iteration 1: left -> right
iteration 2: right -> left
iteration 3: left -> right
所以 FR 使用 temperature 限制每一輪的最大步長：
force_norm = sqrt(fx * fx + fy * fy)
step = min(force_norm, temperature)
pos[i] += force[i] / force_norm * step
每輪結束再降溫：
temperature *= cooling_factor
例如：
cooling_factor = 0.95
前幾輪允許大幅整理結構，後幾輪則逐漸縮小步伐：
early iterations
-> move far
-> untangle rough structure
late iterations
-> move carefully
-> reduce oscillation
Temperature 並不保證演算法找到全域最佳解。它只是避免更新太激進，讓 iterative process 比較穩定。
把三部分組起來，一輪 FR 可以寫成：
def fr_iteration(pos_x, pos_y, edges, k, temperature):
num_nodes = len(pos_x)
force_x = [0.0] * num_nodes
force_y = [0.0] * num_nodes
# 1. All-pairs repulsion
for i in range(num_nodes):
for j in range(num_nodes):
if i == j:
continue
add_repulsion(i, j, pos_x, pos_y, k, force_x, force_y)
# 2. Edge-based attraction
for u, v in edges:
add_attraction(u, v, pos_x, pos_y, k, force_x, force_y)
# 3. Apply one temperature-limited step
next_x = list(pos_x)
next_y = list(pos_y)
for i in range(num_nodes):
next_x[i], next_y[i] = move_with_limit(
pos_x[i], pos_y[i],
force_x[i], force_y[i],
temperature,
)
return next_x, next_y
跑 T 次 iterations，時間複雜度大致是：
T * (N * N + E)
在一般 sparse graph 中，E 通常遠小於 N * N，因此 all-pairs repulsion 會主導成本。
case_4 裡有一個只用 Python standard library 的 CPU 範例：
cd case_4
python examples/layout_naive.py \
--nodes 48 \
--edges 96 \
--iterations 50 \
--seed 42
它會建立一張固定 seed 的 connected graph，從隨機 positions 開始，執行 50 輪同步 FR iterations，再用相同的 proper-crossing 定義比較前後結果。
這一組輸出是：
before
K = 55
crossing pairs = 1,141
mean edge length = 524.49
min node distance = 17.60
after
K = 15
crossing pairs = 197
mean edge length = 274.01
min node distance = 54.02
第一輪最大 displacement 被 temperature 限制在 100.0，第 50 輪降到約 8.10。這讓我們同時看到 attraction、repulsion 與 cooling 的效果：相連 nodes 平均更靠近，nodes 之間的最小距離反而拉大。
這次 K 也從 55 降到 15，但這只是單一 graph 與 seed 的結果。範例仍在 continuous coordinates 上計算，沒有經過 integer quantization、V1～V4 repair 或 official admission，所以不能把 K = 15 當成可提交成績。
這個 baseline 的用途是固定演算法語意，讓 Day 25 的 GPU kernel 能和同一組 force update 對照。
上面的程式先算完所有 forces，再把結果寫到 next_x、next_y。
如果直接在計算途中修改 pos[i]：
for i in range(num_nodes):
compute_force(i)
pos[i] += displacement[i]
後面的 node 會讀到一半新、一半舊的 positions：
node 0 force -> based on old positions
update node 0
node 1 force -> sees new node 0, but old node 2 ... N-1
結果會依賴 node iteration order，也很難和 GPU 平行版本對照。
正確的同步 iteration 應該是：
current positions
|
+-> compute all repulsion
+-> compute all attraction
|
v
displacement buffer
|
v
next positions
GPU kernel 天然會讓很多 threads 同時工作，但這不代表 race condition 自動消失。計算 force 的 kernels 必須只讀 current positions，等它們完成後，另一個 phase 才能更新 positions。
先不做 shared-memory tiling，最直覺的 repulsion kernel 是一個 thread 負責一個 node：
__global__ void repulsive_forces(
const float* pos_x,
const float* pos_y,
float* force_x,
float* force_y,
int num_nodes,
float k_sq
) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i >= num_nodes) return;
float px = pos_x[i];
float py = pos_y[i];
float fx = 0.0f;
float fy = 0.0f;
for (int j = 0; j < num_nodes; ++j) {
if (i == j) continue;
float dx = px - pos_x[j];
float dy = py - pos_y[j];
float dist_sq = dx * dx + dy * dy;
if (dist_sq > 1e-8f) {
fx += dx * k_sq / dist_sq;
fy += dy * k_sq / dist_sq;
}
}
force_x[i] = fx;
force_y[i] = fy;
}
它和 Day 22 的 full crossing kernel 有點像：
one thread owns one output node
one local accumulator in registers
one final global-memory write
Attraction 則適合一個 thread 處理一條 edge：
thread edge_id
-> read edge_u[edge_id], edge_v[edge_id]
-> calculate spring force
-> update force[u]
-> update force[v]
但許多 edges 可能共享同一個 node，所以需要：
atomicAdd(&force_x[u], fx);
atomicAdd(&force_y[u], fy);
atomicAdd(&force_x[v], -fx);
atomicAdd(&force_y[v], -fy);
最後再用一個 thread 更新一個 node：
repulsion kernel
|
attraction kernel
|
v
force buffers
|
v
apply-displacement kernel
|
v
next positions
這已經有 GPU parallelism，但還不是深度優化版。
Day 22 的 crossing kernel 為了避免熱門 edge counter 競爭，選擇重複計算 (e, f) 與 (f, e)。
Attraction 也可以讓每個 node 掃自己的 adjacency，完全避開 atomics：
thread node u
-> scan every incident edge
-> accumulate only force[u]
但 edge-based mapping 也有優點：每條 edge 只計算一次 distance 和 force，然後把相反方向的結果加到兩個 endpoints。
edge-based
less duplicated arithmetic
needs atomic scatter
node-based
single-writer output
edge contribution may be calculated twice
哪一個比較快，取決於 degree distribution。
如果少數 hub nodes 接了大量 edges，atomic contention 可能很嚴重；如果 graph degree 低而且分散，edge-based atomics 可能已經足夠便宜。
所以資深工程師不會只看到 atomicAdd 就立刻重寫。他會先量 degree、contention 與 attraction 在整個 iteration 中的時間占比。因為更大的瓶頸通常還在 N * N repulsion。
Repulsion kernel 中，每個 thread 都會掃過：
pos[0]
pos[1]
pos[2]
...
pos[N - 1]
同一個 block 裡的 threads 其實需要相同的一批 partner positions：
thread i=0 reads partner positions 0 ... N-1
thread i=1 reads partner positions 0 ... N-1
thread i=2 reads partner positions 0 ... N-1
...
但 naive kernel 每個 thread 都從 global-memory view 重新讀取。
GPU cache 可能幫忙，所以「程式碼寫了 N 次 load」不等於一定真的打到 DRAM N 次；不過這個 reuse pattern 非常明顯，值得在下一版主動放進 shared memory。
一批 partner positions
-> block cooperative load once
-> every thread reuses the tile
這會是 Day 25 的主題。
N * N * 2 Tensor在 Python / CuPy 裡，all-pairs force 很容易寫成 broadcasting：
diff = pos[:, None, :] - pos[None, :, :]
程式只有一行，看起來很漂亮，但 diff 的 shape 是：
(N, N, 2)
如果使用 float32，光是這個 tensor 就需要：
memory = N * N * 2 * 4 bytes
當 N = 100,000：
100,000 * 100,000 * 2 * 4
= 80,000,000,000 bytes
= 80 GB
而這還沒有算 distance、force magnitude、direction 與其他 temporaries。
專案的 CuPy 版本會在大圖使用 chunking，避免一次 materialize 完整 N * N；CUDA backend 則讓 thread 直接累加 contribution，只保留 O(N) 的 position 與 force buffers。
這也是一個很實用的 GPU 經驗：高階 array expression 很短，不代表它的 memory footprint 很小。
Force iterations 通常在連續座標中工作：
node 7 -> (0.18374, 0.92751)
但 LCN instance 使用整數畫布，而且 Day 21 定義了 V1～V4：
V1: node 在畫布內
V2: nodes 不共用座標
V3: node 不落在其他 edge 內部
V4: edges 不重疊成一段
因此 layout initializer 的輸出 pipeline 是：
float force positions
|
v
rescale to canvas
|
v
round to integer grid
|
v
resolve duplicate positions
|
v
repair V3 / V4 if needed
|
v
official validation
|
v
exact crossing score
專案的 quantize_to_grid(...) 會做線性縮放、round，以及用 expanding-square search 解決重複座標。這能保證 bounds 和 unique positions，卻不能單獨保證 V3、V4。
所以不能拿 float force energy 當作最終成績。最後仍然要 quantize、repair、validate，再用 Crossing evaluator 算真正的 K。
Force-directed layout 最小化的不是 K。
它偏好的是：
nodes 不要太近
connected nodes 不要太遠
整體幾何分布比較均勻
Crossing 只是這些 forces 可能間接改善的現象，沒有直接出現在 FR 公式裡。
原專案留有一份很小的 layout-init A/B 紀錄。它只有兩個 instances、單一 seed，因此不能拿來做一般性的演算法排名；但很適合提醒我們不要把 force convergence 和 LCN quality 當成同一件事。
| Instance | Strategy | Init time | Initial K | Final K after 15s SA |
|---|---|---|---|---|
| instance_01, 150 nodes | spectral | 3.440 s | 159 | 0 |
| instance_01, 150 nodes | fmmm_lite | 0.480 s | 813 | 0 |
| instance_05, 500 nodes | spectral | 0.479 s | 155,081 | 276 |
| instance_05, 500 nodes | fmmm_lite | 0.110 s | 262,438 | 297 |
在這兩筆資料裡，FMMM-lite 初始化比較快，initial K 卻比較差；經過相同 15 秒 SA 後，instance_01 都到 K = 0，instance_05 則仍由 spectral seed 稍微領先。
這些數據只說明這兩個 instances 和這一個 seed 的結果。它們不能證明 spectral 永遠比較好，也不能證明 force layout 沒有價值。
它真正告訴我們的是：
layout initializer 的速度
layout 看起來是否整齊
initial K
固定 SA 預算後的 final K
是四個不同的測量問題
一個 initializer 也可能 initial K 普通，卻把 graph 放進比較容易搜尋的 basin，最後反而得到更好的結果。因此正式比較應該使用多個 instances、seeds，並固定 repair 與 SA budget。
Crossing 和 Layout 都有 all-pairs 工作，看起來很像：
Crossing
edge pairs
Layout repulsion
node pairs
但它們的最佳化方向不同。
Crossing 的關鍵是 moved node 只影響少量 incident edges：
先從演算法上刪掉不可能改變的 pairs
Layout repulsion 每一對 nodes 都真的有貢獻，不能直接說遠處 node 一定不用算：
先讓同一批 positions 在 GPU 上被有效重用
這個案例會帶出另一組 GPU 觀念：
global memory traffic
shared-memory tiling
register accumulation
coalesced load
kernel fusion
small-graph 與 large-graph backend
也會帶出一個和 Day 23 不同的限制：演算法輸出的是 initializer，不是 objective 的 exact answer。算得快之外，還要追蹤 downstream 的 validated K。
Day 24 完成的是 force-directed layout 的基本模型：
N * N，通常主導 iteration。K 沒有單調關係，必須用 downstream 結果評估。今天的 kernel 已經能讓很多 nodes 同時計算，但每個 block 仍然反覆讀取同一批 partner positions。
Day 25 要解的就是這筆記憶體帳：
同一個 position tile
能不能只從 global memory 載入一次，
再讓整個 block 的 threads 重用？
答案會把我們帶到 shared-memory tiling、每個 thread 擁有多個 nodes、迭代融合，以及一個反直覺結果：融合越多，register 壓力越高，occupancy 不一定越好。
(等我放了我再更新連結 現在還沒有檔案XD 26-09-13)
iThome鐵人賽
看影片追技術